Universal slow growth of entanglement in interacting strongly disordered systems 
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Recent numerical work by Bardarson et. al. [Phys. Rev. Lett. 109, 017202 (2012)] revealed a 
slow, logarithmic in time, growth of entanglement entropy for initial product states in a putative 
many-body localized phase. We show that this surprising phenomenon results from the dephasing 
due to exponentially small interaction-induced corrections to the eigenenergies of different states. 
For weak interactions, we find that the entanglement entropy grows as ^\ii(Vt/H), where V is the 
interaction strength, and £ is the single-particle localization length. The saturated value of the 
entanglement entropy at long times is determined by the participation ratios of the initial state over 
the eigenstates of the subsystem. The proposed mechanism is illustrated with numerical simula- 
tions of small systems. Our work shows that the logarithmic entanglement growth is a universal 
phenomenon characteristic of the many-body localized phase in any number of spatial dimensions, 
and reveals a broad hierarchy of dephasing time scales present in such a phase. 

PACS numbers: 72.15.Rn, 05.30.Rt, 37.10.Jk 



VH 

a 

a 
o 
o t 

<N 
> 

m 

o 
\o 
ti- 

o 
en 



s 



Introduction. While it is well-known that arbitrar- 
ily weak disorder localizes all single-particle quantum- 
mechanical states in 1 and 2 dimensions, the effect of 
disorder potential on the states of interacting many-body 
systems largely remains an open problem. Works [2, 3] 
conjectured that localization in a system of interacting 
electrons survives in the presence of weak interactions. 
When the strength of the interactions is increased, at 
some critical value a transition to the delocalized phase 
("many-body localization transition") takes place. The 
existence of such a transition was further supported by 
numerical simulations of small systems [4-8, 12-16]. 

An important challenge is to understand the physi- 
cal properties of the localized phase of interacting parti- 
cles. Recent numerical work [1] (also see Ref. [5]) revealed 
that even very weak interactions dramatically change the 
growth of entanglement of non-equilibrium many-body 
states. Authors of Ref. [1] studied the time evolution of 
product states in a one-dimensional disordered XXZ spin 
chain, which can be reformulated in terms of interacting 
fermions in a random potential (the same model was also 
considered by other authors in the context of the many- 
body localization [8]). In the absence of interactions, 
such states maintain a low degree of entanglement upon 
evolution, and the entanglement entropy obeys an area 
law. In contrast, in the presence of interactions the states 
showed a slow, logarithmic in time, growth of entangle- 
ment entropy (here and below we use "entanglement" and 
"entanglement entropy" interchangeably). The saturated 
value of the entanglement entropy was found to vary ap- 
proximately linearly with system size, and remained well 
below the maximum possible entanglement entropy of the 
system [1, 9] 

In this paper, we identify a mechanism that underlies 
the logarithmic growth of entanglement in interacting lo- 



calized states. The key observation is that although very 
weak interactions affect the many-body eigenstates of a 
non-interacting Hamiltonian only weakly, they do change 
the energies of these states. It is these small interaction- 
induced corrections to the energy which lead to the de- 
phasing between different eigenstates on long time scales. 
We argue that this gives rise to a logarithmic growth of 
entanglement with time for a broad class of initial states 
that are a product of states in the two subsystems, in- 
cluding, but not limited to the initial states considered 
in Ref. [1]. 

For weak interactions, our mechanism leads to the fol- 
lowing predictions regarding entanglement growth as a 
function of system's parameters: (i) entanglement en- 
tropy grows as S en t(t) oc £\og(Vt/h), where V is the in- 
teraction strength and £ is the single-particle localization 
length; (ii) the saturation value of S ent is of the order of 
"diagonal entropy" Sdiag [17] of a subsystem for a given 
initial state. Diagonal entropy is determined by the par- 
ticipation ratios of initial states in the basis of single- 
particle eigenstates. We also illustrate these predictions 
with numerical simulations of small systems, in particular 
by constructing examples of initial states for which the 
saturated entanglement is equal to the diagonal entropy. 

Model. Without loss of generality, we consider a 
ID lattice model of fermions with on-site disorder and 
nearest-neighbor interactions, 



H = t 2_. a l a j + /_^ Wialai + V Y^ fiihj, 



,t. 



(ij) 



(1) 

where i,j = 1, ...N, and (ij) denotes nearest neighbors. 
This model is equivalent to the XXZ spin chain with a 
random field along the z-direction considered in Ref. [1]. 
From our discussion below, it will become clear that the 
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FIG. 1. (Color online) (a) Generation of entanglement for the 
case of a single particle in each half of the system prepared 
in a superposition of two eigenstates. Exponentially small 
overlap of the orbitals leads to the dephasing time growing 
exponentially with distance between particles, (b) Entangle- 
ment entropy as a function of time for a given realization of 
disorder and different interaction strengths. When V = 0, 
S e nt ~ 10 -4 and remains small for all times (for the plot 
tmax = 2.5 • 10 11 ). For non-zero interactions, S e nt(t) collapse 
on a single curve when time is scaled by 1/V. System size is 
L — 10 sites, disorder strength is W — 6. Initially, particles 
are in equal superposition of two sites at the opposite edges 
of the system [see panel (a)]. As these states are not the ex- 
act eigenstates, the maximum value of S eVL t is slightly below 
ln2^0.69. 



logarithmic growth of entanglement in interacting local- 
ized systems is a robust phenomenon which does not de- 
pend on the microscopic details, nor on the dimensional- 
ity of the system. 

Throughout the paper, we focus mostly on the regime 
of weak interactions, for which the logarithmic growth of 
entanglement found in [1] is perhaps the most striking. 
In the absence of interactions, V = 0, disorder local- 
izes the single-particle states, with localization length £. 
The many-body eigenstates are simply states in which a 
certain number of single-particle orbitals is occupied. In- 
teractions that are much weaker compared to the typical 
level spacing on the scale £ do not significantly modify 
the many-body eigenstates. We have explicitly verified 
this statement for small systems, and assume it holds in 
general. However, even though the eigenstates are not 
strongly affected by the interactions, their energies are 
modified. Interactions between remote particles yield a 
contribution to the energy that is exponentially small in 
the ratio x/£, where x is the distance between particles. 
Effectively, if we fix the positions of all particles, except 
for a pair of particles which are situated at a distance 
x ^> £ away from each other, the interaction energy of 
this pair is of the order Ve~ x ^^ and the corresponding 
dephasing time is tdeph ~ he x >^ /V. This gives rise to a 
hierarchy of dephasing time scales present in the prob- 



lem, ranging from the fastest t m [ n = H/V to the slowest 
t max = t m i n e L /^, where L is the system size. 

Generally, the product initial states considered in 
Ref. [1], as well as the initial states of other kinds con- 
sidered below, are a superposition of many eigenstates. 
The interactions introduce a slow dephasing between dif- 
ferent states, and effectively generate entanglement be- 
tween different remote parts of the system. A subsys- 
tem of size x becomes nearly maximally entangled with 
the rest of the system after an exponentially long time 
^deph(^) ~ he x '^ /V . Thus, if one considers a bipartite 
division of the system, the entanglement entropy will in- 
crease logarithmically in time, as more and more degrees 
of freedom in the two subsystems become entangled with 
each other. 

Two particles. Let us start with a simple example 
which demonstrates that the slow growth of entangle- 
ment occurs for just two particles. More specifically, con- 
sider two distant particles prepared in an equal- weight su- 
perposition of two neighboring localized orbitals, |^o) — 

Cl /^ 2 ° 3 J; 4 |Q) ; where c\ creates an eigenstate localized 
near site i. We assume that the distance between wave 
functions 1,2 and 3,4 is large, x ^> £. This situation is 
depicted in Fig. 1. 

In the absence of interactions, no entanglement would 
be generated during time evolution. Interactions change 
this picture by introducing a correction to the energy of 
the state \a/3) = c^cUO), where a = 1, 2, /? = 3,4. In the 
leading order of perturbation theory, the energy of this 
state is given by 

E*p =Ea+ep + 5E a(3 , 5E a(3 = C a(3 Ve- x/ ^ (2) 

where e a , ep are the single-particle energies, the last term 
is introduced by the interactions, and C a p is a constant 
which depends only algebraically on x. 
The time-evolved state is given by 

a=l£=3 

and the reduced density matrix for the first particle reads 



1 



PL 



1 



F(t)/2 \ 
2 \F*(t)/2 * J ' 



(4) 



where F(t) = e~ int (l + e~ ibm ), SQ = SE 14 - 5E 2A - 
5Ei3 + 5E 23 , and Q = e\ - e 2 + SEis - SE 23 . The eigen- 
states of pl therefore oscillate with a very long period 
T = 27T/6Q ~ (h/V)e x ^. At times t = (2n + 1)tt/£Q, 
the off-diagonal matrix elements vanish, and eigenstates 
become equal to 1/2. At these times, the two particles 
become maximally entangled with S ent = In 2. Fig. 1 
demonstrates that even weak interactions lead to the en- 
tanglement of the order of 5 en t ~ In 2, and the rate of 
entanglement change is inversely proportional to the in- 
teraction strength. Note that no disorder or time aver- 
aging is used. 



General case. Turning to the general many-body 
case, let us divide the system into two parts, C and 7Z, la- 
beling the single-particle orbitals that are localized dom- 
inantly in the left part by index a n , and those residing 
in the right part by j3 n . There may be some ambiguity 
for the state residing near the boundary between C and 
7£, but we will be interested in systems of size L ^> £, for 
which the boundary effects are not very important. 

We consider initial states that are products of some 
superposition of states with a definite number of particles 
in C and 1Z: 



l*(* = 0))= Y, A {a}\ai-aK) 

{a}e£ 



^ B m \^.^ M ). 
Wen 

(5) 
Coefficients A, B are chosen such that \I/ is normalized. 

Neglecting the change to the eigenstate due to inter- 
actions, the reduced density matrix for C after time evo- 
lution reads 



pc = y^Paa>\a)(a'l 



(6) 



Paa> = A a A* a , J2 \ B ?\ 2 ex P [i( E a'j3 ~ E a/3 )t] , (7) 

where we used a short-hand notation a = {a},/3 = {/?}. 
It is convenient to define A a (t) = A a e~ lEat , where E a 
is the energy of \a) state for the isolated C subsystem. 
Assuming that \a) x |/3) remains an eigenstate (this may 
not be true near the boundary, but the boundary effect 
is not important for entanglement growth, at least in 
large systems), the above equation, written in terms of 
coefficients A a (t), preserves the same form, except the 
energies E a p should be substituted by the interaction 
energy SE a p between particles in C and 1Z subsystems. 
For particles that reside far away from the boundary, this 
correction can be calculated in perturbation theory. 

The energy difference 5E a >p — 5E a p that enters the off- 
diagonal elements of the reduced density matrix, to the 
leading order, is proportional to Ve~ x ^. Here x is the 
minimum distance between a particle in £, the position 
of which is different in states a and a', and the particles 
in 1Z. However, it also contains many smaller contri- 
butions, which arise due to the interaction between more 
distant pairs of particles. Thus, the off-diagonal elements 
oscillate at a number of very different, incommensurate 
frequencies. 

The interaction energy leads to dephasing, which de- 
creases the off-diagonal elements of the density ma- 
trix, thus generating entanglement. Effectively, at times 
t(x) ~ toe x ^ the degrees of freedom within a distance 
x(t) ~ log (t /to) from the boundary between C and 1Z 
are affected by dephasing, while states that differ only in 
the positions of particles further away from the boundary 
are still phase-coherent. This generates the entanglement 
entropy given by 

Sent(t) = CS diag , S diag = -J2P i (x)lnP i (x), (8) 



where Pi{x) are the probabilities of different states \a) 
in a segment of size x, calculated from the wave func- 
tion of the initial state. Quantity Sdiag is the diagonal 
entropy - a maximum achievable entropy for a given ini- 
tial state, assuming that interactions do not change the 
eigenstates. S en t is expected to be smaller than Sdiag 
by a factor C < 1; the precise value of this pre- factor is 
non-universal, and depends on the preparation of the ini- 
tial state. In the long-time limit, assuming that 1Z ^> C 
and the initial state is a superposition of many differ- 
ent states, we expect the off-diagonal elements to be- 
come very small such that the entanglement entropy ap- 
proaches its maximum value with C — > 1. 

Since for initial product states Sdiag is proportional to 
the subsystem size, entanglement grows logarithmically: 



Sent(t)oc£log(Vtyft). 



(9) 



We emphasize that our argument does not rely on av- 
eraging, and therefore entanglement grows according to 
Eq.(9) even for a single disorder realization, and even for 
relatively small systems. 

Numerical simulations. In order to illustrate the 
above mechanism, and to explore the growth of entangle- 
ment for different initial states, we performed numerical 
simulations of finite systems. We use XXZ spin chain 
as an effective description of fermions in ID, with a ran- 
dom magnetic field uniformly distributed in the inter- 
val [— W, W]. Hopping, described by exchange J_l, is set 
to one, whereas the interaction strength is controlled by 
J z = V . In the many-particle simulations, the Hilbert 
space is restricted to the total S z = sector. Using exact 
diagonalization for systems up to 14 spins, we compute 
the time evolution of various initial states of the system. 
This allows us to obtain the entanglement entropy S^nt (t) 
between two halves of the system as a function of time 
(even number of sites is considered). We study its aver- 
age S ent (t) over different initial states belonging to the 
same class, and over different realizations of disorder. 

First, similar to Ref. [1], we consider a class of prod- 
uct states of up or down spin on every site ( "up-down" 
states). Results for S ent (t) for a system of L = 12 sites 
and disorder W = 5 are shown in Fig. 2(a). After a rapid 
increase of entropy during time t ~ 1/Jj_ due to diffu- 
sive transport on a scale smaller than the localization 
length, S en t(t) saturates for a non-interacting system. In 
the presence of even weak interactions, S ent (t) continues 
to grow further. The growth rate is found to be propor- 
tional to ]n(Vt/h) [Fig. 2(b)], in full agreement with our 
analysis above. 

The saturation value of the entanglement, *S en t(oo), 
does not vary appreciably with the interaction strength 
when interactions are weak. This further supports the 
conclusion that weak interactions only weakly alter the 
eigenstates of the system. For fixed V = 0.01, S ent (oo) 
decreases with disorder W [Fig. 2(c)]: for stronger disor- 
der, the initial up-down state becomes increasingly more 
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FIG. 2. (Color online) (a) Averaged entanglement entropy of 
a random "up-down" state as a function of time shows a char- 
acteristic logarithmic growth with time on long time scales, 
(b) S en t — So collapses onto a single curve as a function of 
hi(Vt/h). So is the saturation entropy for a non-interacting 
system. System size L = 12, W = 5. (c) Saturated en- 
tanglement for V — 0.01 decreases as a function of disorder 
strength, (d) Ratio C = Sent (oo)/ Sdiag shows a slower de- 
crease with W. Error bars (if not shown) are approximately 
equal to the size of symbols in each plot, with the exception 
of panel (b) , where they are omitted for clarity. 
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FIG. 3. (Color online) Saturation entropy as a function of 
disorder W for a strongly entangled state (a), and a product 
of a strongly entangled state and an "up-down" state (c). 
In the former case, the saturated entropy is approximately 
independent of W, while in the latter it decreases as W is 
increased. The ratio between saturation and diagonal entropy 
for two kinds of initial states, (b) and (d). In the second 
case, this ratio tends to one for larger system sizes as W is 
increased. Interaction strength is set to V = 0.01. Error bars 
(if not shown) are approximately equal to the size of symbols 
in each plot. 



similar to an eigenstate of the system, and its diagonal 
entropy diminishes. 

We compare the saturated entanglement to Sdi ag , cal- 
culated using the values of Pi(L/2) for the initial state 
of the C subsystem. The Pi(L/2) are obtained from the 
density matrix, using the eigenstates of the interacting 
Hamiltonian restricted to C. In this sense, while S en t(oo) 
is determined from the time evolution of the system, Sdiag 
is solely the property of the initial state. 

To interpret the dependence of the ratio C = 
S en t(oo) / Sdiag on system size and disorder, Fig. 2(d), we 
must take into account two additional effects important 
for small systems: (i) the diffusion of particles across the 
entanglement cut; and (ii) the inefficiency of decoherence 
when the number of terms in Eq. (5) is small or when C 
and 1Z are of equal size. These effects counteract each 
other, as diffusion leads to an additional contribution 
to the entanglement not captured by (8). On the other 
hand, inefficient decoherence leads to incomplete dephas- 
ing, and decreases the saturated entanglement compared 
to Sdiag- The positive contribution from (i) is suppressed 
for larger systems or smaller localization lengths. The 
effect of (ii) depends on the initial state. For the "up- 



down" states in the localized phase, the participation ra- 
tio is on order of unity, and the effect of (ii) is very pro- 
nounced. Thus, C is smaller than one, and it decreases 
with stronger disorder or system size [Fig. 2(d)]. 

Next, we consider another kind of initial states, which 
have larger participation ratios. The initial state of 1Z 
and C subsystem is chosen as a projection to S z = 
sector of the product state of (| t) ± | \))/y/2 per each 
site, with ± signs chosen at random. Spins within each 
subsystem are therefore strongly entangled, but there is 
no entanglement between the subsystems at t = 0. We 
found that the logarithmic entanglement growth occurs 
also in this case, but the saturated entropy, Fig. 3(a), is 
larger compared to the previous case, and varies weakly 
with disorder. Ratio of saturated to diagonal entropy, 
Fig. 3(b), now scales to one when the system size is in- 
creased, contrary to the "up-down" initial states. For 
this type of initial state, due to larger values of 5 en t(oo), 
boundary diffusion contribution to entanglement is less 
important; also, the superposition of a large number of 
eigenstates in each half of the system makes decoherence 
more efficient, thus r is closer to one. 

Finally, we construct an example where Sent reaches 



^diag- We take a product of the "up-down" state in 
£, and the strongly entangled state in 1Z. To further 
suppress the diffusion, we require the two sites adjacent 
to the entanglement cut to be always in the | \) state. 
'S'ent(oo) displays the behavior similar to the "up-down" 
case, Fig. 3(c), but is larger due to the more effective 
dephasing. Remarkably, Fig. 3(d) demonstrates that for 
larger system sizes, saturation and diagonal entropies be- 
come equal, in agreement with the above analysis. 

Discussion. To summarize, we presented a mecha- 
nism for the logarithmic growth of entanglement in a 
model undergoing a many-body localization transition. 
We also established the laws that govern the entangle- 
ment growth, and tested them in numerical simulations 
for different initial states. 

Although we focused on the limit of weak interactions, 
in which the eigenstates are similar to those of a non- 
interacting model, we expect our conclusions to hold 
also for stronger interactions which locally do modify the 
eigenstates. In this case, the saturated value of the entan- 
glement is expected to be determined by the participation 
ratios of the initial state in the basis of interacting sub- 
system's eigenstates. Furthermore, our conclusions are 
expected to apply to localized interacting systems in any 
number of spatial dimensions 

Our work also points to an exponentially broad distri- 
bution of dephasing time scales present in a many-body 
localized system. This has interesting implications for co- 
herent dynamics of a small quantum system interacting 
with a many-body localized environment. For simplicity, 
let us consider a situation where the small quantum sys- 
tem is simply one of the particles that constitute the sys- 
tem, prepared in a superposition of some single-particle 
orbitals. Then, if a region of size L* surrounding this 
particle is prepared in an eigenstate of the Hamiltonian 
(while the rest of the system is presumably in a mixed 
state), the dephasing time of the particle will be pro- 
portional to exp(L*/£) and therefore exponentially long. 
This suggests that in certain systems it might be possible 
to reduce decoherence by introducing strong disorder in 
the environment. 

Finally, we note that recently Vosk and Alt man [18] 
considered an XXZ model with random exchange inter- 
action, but without random field. For a special initial 
state, they developed a strong-disorder renormalization 
group (RG) procedure, and computed entanglement en- 
tropy, finding that it grows as a power of logt. The dif- 
ference from our result stems from the fact that the state 
considered in Ref. [18] was critical, however, the basic 
underlying mechanism - dephasing due to exponentially 
weak interactions between remote spins - is qualitatively 
similar. 

Prior to submission, we became aware of a related in- 
dependent work by Huse and Oganesyan [19]. Using a 



phenomenological description of the many body local- 
ized phase, they also predicted a logarithmic growth of 
entanglement with time. Our results agree insofar as they 
overlap. 
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